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ABSTRACT 


A  statistical  procedure  for  determining  the  resolution  to  be  used  in  spectral- 
matrix  estimation  is  needed  for  the  exploratory  analysis  of  non-stationary  multiple 
time  series.  Our  procedure  accomplishes  this  by  creating  periodograms  for  a 
frequency-time  interval  and  then  successively  dividing  this  interval  into  subintervals. 
At  each  step,  the  division  is  chosen  by  comparing  spectral  matrices  using  the 
largest-root  statistic.  This  multiple-decision  rule  is  shown  to  be  consistent. 
Small-sample  biases  and  approximate  critical  regions  for  testing  stationarity  and 
frequency  uniformity  are  found.  Simulation  shows  that  simultaneous  analysis  of 
all  the  series  performs  better  than  individual  analyses. 
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1.  INTRODUCTION 


In  all  spectral  analysis,  a  key  problem  is  choosing  the  resolution  so  that  the 
spectral  estimates  have  adequate  stability  without  unacceptable  bias.  Compared 
to  the  case  of  a  single,  stationary  time  series,  this  problem  is  significantly  more 
complicated  in  the  case  of  non-stationary  or  multiple  time  series.  In  particular, 
the  usual  techniques  for  determining  the  proper  frequency  or  time  resolution  from 
the  data  are  not  effective  for  non-stationary  multiple  time  series.  In  this  paper, 
we  analyze  a  technique  suitable  for  non-stationary  multiple  time  series  which  we 
proposed  previously  [20]  . 

The  time  series  for  which  our  procedure  is  intended  are  those  with  slowly- 
varying  covariance  properties  and  those  for  which  spectral  analysis  is  appropriate. 
Slowly  varying  means  that  in  effect  the  covariance  matrix  depends  only  on  the  time 
difference  within  time  intervals  over  which  statistical  dependence  exists.  This 
notion,  which  justifies  the  estimation  of  the  spectral  matrix  as  a  function  of  time, 
has  been  made  more  precise  by  Priestley  [25]  .  The  time  series  for  which  spec¬ 
tral  analysis  is  appropriate  are  generally  large,  have  zero  mean  or  the  mean 
removed,  and  have  characteristics  of  Interest  that  are  most  naturally  expressed 
in  the  frequency  domain.  Such  processes  have  been  considered  in  the  statistics 
literature  in.  25,  28] .  Variation  of  the  spectral  properties  with  time  is  often 
of  interest  in  applications. 

Spectral-matrix  estimates  are  rarely  the  end  result;  they  usually  must  be 
interpreted.  Beside  graphical  interpretation  of  the  spectrum,  multiple  time  series 
can  be  analyzed  by  applying  to  the  spectral  matrix  techniques  useful  for  interpreting 
covariance  matrices  such  as  principal-component  analysis  [2,  3,  10,  12,  27], 

The  bias-stability  trade-off  for  non-stationary  multiple  time  series  differs  in 
three  important  ways  from  the  trade-off  for  a  single  stationary  time  series.  Gener¬ 
ally,  these  differences  lead  to  requirements  for  both  more  resolution  and  more 
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stability  in  the  former  case.  The  first  difference  is  that  for  non-stationary  pro¬ 
cesses,  inadequate  time  resolution  introduces  smearing  in  time  that  is  equivalent 
to  the  smearing  introduced  by  inadequate  frequency  resolution.  If  the  spectral 
estimate  is  formed  from  data  in  a  time  interval  during  which  the  population  spec¬ 
tral  matrix  changes,  the  mean  of  the  estimate  will  be  an  average  of  the  population 
spectral  matrix  over  the  time  interval.  To  understand  the  second  difference,  con¬ 
sider  a  time-varying  linear  filter  with  one  input  and  several  outputs.  Since  there 
is  only  one  input,  the  population  spectral  matrix  of  the  outputs  will  have  rank  one. 
However,  a  spectral-matrix  estimate  obtained  from  an  interval  in  which  the  rela¬ 
tions  among  the  outputs  vary  with  frequency  or  time  will  not  have  rank  one.  This 
is  the  bias  that  can  cause  a  coherence  estimate  that  should  be  one  to  be  much  less 
than  one.  Analogous  smearing  can  occur  with  more  complicated  spectral-matrix 
structures.  The  third  difference  is  the  requirement  of  the  techniques  for  spectral- 
matrix  interpretation  for  more  stability  than  is  needed  for  interpretation  of  the 
spectrum.  The  procedures  of  complex  multivariate  analysis  have  biases  due  to 
inadequate  stability  because  they  form  non-linear  functions  of  the  spectral-matrix 
estimates.  For  example,  the  usual  coherence  estimates  are^biased  upward  due  to 
inadequate  stability.  Thus,  in  coherence  estimation,  there  is  a  trade-off  between 
two  biases  which  act  in  opposite  directions  [30]  . 

The  procedures  that  have  been  proposed  for  determining  the  proper  resolution 
from  the  data  fall  into  four  categories,  those  based  on  visual  interpretation  of  data 
displayed  versus  frequency  and  time,  those  that  test  for  serial  correlation,  those 
that  determine  the  proper  frequency  resolution  for  stationary  time  series,  and 
those  that  test  for  stationarity.  Our  procedure  differs  from  each  of  these.  Although 
we  have  not  done  a  thorough  performance  comparison,  our  procedure  seems  to 
have  unique  capabilities  useful  in  many  situations  including  the  exploratory  analy¬ 
sis  of  multiple  time  series  with  many  components. 

A  time-honored  approach  to  determining  the  frequency  and  time  resolution  for 
a  single  time  series  is  the  spectrogram  [18,  29],  A  spectrogram  is  created  by 
calculating  the  periodograms  for  successive  time  intervals  and  then  displaying 
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them  as  intensity  versus  frequency  and  time.  As  a  tool  for  exploratory  analysis, 
the  spectrogram  is  undoubtedly  effective  for  a  single  time  series.  Extension  to 

multiple  time  series  is  largely  blocked  by  the  inability  to  simultaneously  display 
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the  p  independent  variables  (for  p  time  series)  of  the  spectral  matrix  versus  fre¬ 
quency  and  time. 

Our  procedure  is  a  generalization  of  the  test  for  serial  correlation  in  a  single 
stationary  time  series  presented  by  Durbin  [8].  This  test  compares  the  sample 
path  of  the  cumulated  periodogram  with  Kolmogorov-Smirnov  limits.  Section  3 
clarifies  this  relation. 

For  stationary  time  series,  the  technique  of  trying  different  frequency  reso¬ 
lutions  has  often  been  suggested.  One  procedure  is  a  graphical  comparison  of 
spectra  with  different  resolutions  [15,  p.  280] .  In  autoregressive  spectral  esti¬ 
mation,  a  more  exact  procedure  is  possible  for  choosing  the  proper  order  for  the 
autoregressive  scheme  and  thus  the  resolution  [l,  23].  This  latter  procedure 
extends  to  multiple  time  series.  Except  for  special  types  of  nonstationarity,  these 
procedures  have  not  been  extended  to  finding  the  time  resolution  for  non-stationary 
time  series. 

Tests  for  stationarity  have  been  proposed  [22,  p.  405;  26;  28],  These  tests 
are  based  on  comparison  of  the  results  from  non-overlapping  frequency-time 
intervals  that  are  large  enough  to  allow  stable  spectral-matrix  estimates  to  be 
obtained.  Reducing  the  problem  to  comparing  results  that  are  independent  and 
identically  distributed  under  the  null  hypothesis  has  great  advantages  including 
simplicity  and  robustness.  For  example,  Priestley  reduces  the  problem  to  the 
analysis  of  variance  [26,  28] .  Unfortunately,  determining  sufficiently  large 
frequency-time  intervals  a  priori  is  not  always  possible.  When  the  spectral  pro¬ 
perties  vary  in  both  frequency  and  time,  intervals  with  the  same  shape  that  are 
sufficiently  large  and  have  constant  spectral  properties  may  not  even  exist.  Thus, 
our  procedure  which  is  capable  of  choosing  intervals  with  a  wider  variety  of  shapes 
is  sometimes  necessary  despite  the  increased  complexity.  This  is  especially  true 
in  exploratory  studies  where  uncovering  the  unexpected  is  desired  [7,  3l]. 
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In  the  following  sections,  we  describe  our  procedure,  discuss  its  statistical 
properties,  and  report  the  results  of  simulation.  In  Section  2,  we  describe  the 
rationale  for  our  procedure  and  some  alternatives  which  might  be  useful.  Analy¬ 
sis  of  the  dependence  of  one  step  of  our  procedure  on  the  previous  ones  does  not 
seem  possible.  Thus,  in  Section  3,  we  derive  statistical  properties  of  one  step 
in  isolation  including  asymptotic  consistency,  univariate  and  approximate  joint 
distributions  for  our  measure  of  dissimilarity  for  the  case  of  constant  spectral 
properties,  and  small-sample  biases.  From  these  results,  approximate  critical 
regions  for  testing  stationarlty  and  uniformity  in  frequency  can  bo  computed.  In 
Section  4,  we  demonstrate  by  simulation  that  analyzing  nn  entire  multiple  time 
series  simultaneously  performs  better  than  analyzing  each  component  time  series 
individually. 
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2.  CLUSTERING  SPECTRAL  ESTIMATES 


This  paper  presents  a  procedure  for  clustering  the  spectral-matrix  estimates 
corresponding  to  different  frequencies  and  times.  By  Fourier  transforming  the 
data  in  successive  time  intervals,  matrix  analogues  of  (modified)  periodogram 
coefficients  can  be  obtained  for  a  grid  of  points  in  frequency  and  time.  These  ele¬ 
mentary  spectral-matrix  estimates  are  approximately  Independent  but  not  identi¬ 
cally  distributed  when  the  population  spectral  matrix  varies  with  frequency  or  time. 
Since  these  estimates  are  unstable,  they  must  be  averaged  to  make  them  useful. 

If  the  elementary  spectral-matrix  estimates  are  nearly  unbiased,  the  problem  of 
obtaining  good  spectral  estimates  is  reduced  to  choosing  frequency-time  intervals 
in  which  the  population  spectral  matrix  is  nearly  constant  and  then  smoothing  over 
these  intervals.  Our  procedure  chooses  rectangular  intervals  by  successive  divi¬ 
sion  of  a  large  rectangular  interval. 

Although  Fourier  transforming  successive  non-overlapping  time  intervals  to 
obtain  periodograms  is  computationally  simple  and  also  the  case  usually  treated  in 
justifying  the  complex-Wishart  distribution,  the  technique  is  known  to  have  rela¬ 
tively  poor  spectral  windows.  An  alternative  with  better  statistical  properties  can 
be  obtained  using  faders  and  overlapping  time  intervals  [4;  12,  p.  265] .  Since  a 
technique  that  is  appropriate  to  our  procedure  has  not  been  detailed  previously,  we 
describe  one  option.  We  assume  that  the  multiple  time  series  has  zero  mean. 

Let  the  sample  of  the  (discrete)  multiple  time  series  be  x(m),  m  =  1, . . . ,  M, 
where  m  indexes  time  and  x(m)  is  the  column  vector  of  the  observations  from  the 
p  component  time  series  at  time  m.  Let  M  =  (T  +  l/2)L  where  T  and  L  are  inte¬ 
gers  so  that  (2T  -  2)  intervals  with  length  2L  and  75%  overlap  can  be  created. 

There  are  three  steps  in  computing  the  elementary  spectral-matrix  estimates. 

First,  from  the  data  in  each  interval  the  finite  Fourier  transforms 


6 


-1/2  O  2L 

w(f',  t»)  =  (4  n  L)  >  ,  x(m  +  t*  L/2)  exp  (-i  n  f'  m/L), 

Lmd  m  =  1 


(2.1) 


1  <  f'  <  L  -  1,  0  <  t'  <  2T  -  3, 


are  computed.  We  exclude  the  exceptional  frequencies  0  and  II.  Second,  a 
cosine-bell  fader  is  applied  to  give 


w(f\  t')  =  -  w(f*  -  1,  t')/4  +  w(f\  t')/2  -  w(f *  +  1,  t')/4 
-1/2  2L 

=  (16nL)  >  x(m  +  t'  L/2)  (1  -  cos  (II  m/L))  (2.2) 

m  =  1 

•  exp(-i  H  f'  m/L),  2  <  f*  <  L  -  2,  0  <  t'  <  2T  -  3. 

Finally,  the  estimates  from  adjacent  frequencies  and  times  are  averaged  to  pro¬ 
duce  the  elementary  spectral-matrix  estimates. 


,  V""*  2f  +  1  2t  +  1  A  A 

a(f,  t)  =  (2/3)2^  f,  =  2f  2-f  f  =  2t  *(f’’  V)  ™*{V'  V)t 


(2.3) 


1  <  f  <  (L  -  3)/2,  0  <  t  <  T  -  2, 


where  *  denotes  conjugate  transpose.  a(f,  t)  is  properly  scaled  to  be  an  estimate  of 


-i  w.  m 


F( 


W_)  =(1/2H)  V*  03  R(m)e  f  ,  w  =  II(4f  +  1)/2L, 
i  m  =  -®  i 


(2.4) 


where  R(m)  is  the  covariance  matrix. 


Comparison  of  the  frequency  window  with  the  time  window  for  our  elementary 
spectral-matrix  estimates  a(f,  t)  shows  considerable  similarity.  The  spectral 
window  of  a(f,  t)  is  given  for  large  L  by 


2  L  cos2  ^L(u)  -  u^)  L^(w  -  w^)2  +  9  n2/4  J 
f  f  3[(L2(w-  Wf)2-9II2/4)(L2(u.-Wf)2-n2/4)J2 


"  Wr)  = 


,  o  <  Wj  <  n. 


(2.5) 
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The  mean  of  a(f,  t)  when  a(f,  t)  is  formed  from  a  stationary  time  series  with  spec¬ 
trum  F(w)  is 

n 

kjw  -  uj J  F(w)  du,  0  <  w  <  II.  (2. 6) 

n  f  r 

For  non -stationary  time  series,  an  analogous  temporal  window  for  a(f,  t)  is  given 
by 

(12  n  L)"1  [l  -  cos  (n  (  (m  -  r)/L  +  5/4))]  2 
for  -  5/4  <  (m  -  t)/L  <  -  3/4 

(12  n  L)"1  j[l  -  cos  (n  (  (m  -  t)/L  +  5/4))]  2 
+  [l  -  cos  (II  (  (m  -  t  )/L  +  3/4))] 
for  -  3/4  <  (m  -  t)/L  <  3/4 

(12  H  L)‘ 1  [l  -  cos  (n  (  (m  -  t)/L  +  3/4))]  2 
for  3/4  <  (m  -  t)/L  ^  5/4 

0  otherwise. 

The  mean  of  a(f,  t)  when  a(f,  t)  is  formed  from  an  uncorrelated  series  with  time- 
varying  variances  R  (0)  is 

v 

E®  k  (m  -r)R  (0),  T  =  t  +  L/4 .  (2. 8) 

m  — -®  t  t  m  t 

The  windows  k  and  k  are  shown  in  Figure  1  along  with  the  windows  for  averages 

X  l 

of  three  adjacent  elementary  spectral  matrix  estimates.  The  sidelobes  on  the  fre¬ 
quency  window  are  very  small.  Within  the  averaging  interval,  the  windows  are 
nearly  flat. 

Our  procedure  is  based  on  the  premise  that  the  a(f,  t)  are  independent.  This 
is  only  approximately  true  for  adjacent  estimates  as  can  be  seen  from  the  moments 
of  a(f,  t)  for  unit-power  white  noise  (RJt)  =  6(t)).  For  large  L,  we  have 


l  (2.7) 
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E  a(f,  t)  =  (2  II) 

var  a(f,  t)  =  (2  II)"2  (0. 527) 

2  (2<9) 
cov  (a(f,  t),  a(f  +  1,  t))  =  (2  II)  (0. 095) 

cov  (a(f,  t),  a(f,  t  +  1))  =  (2  II)"2  (0. 095) . 

If  a  spectral  estimate  is  formed  by  averaging  a(f,  t)  over  f  <  f  <f  ,  t  <  t  <  t 

LULU 

where  f  -  f  and  t  -  t  are  both  large,  then  the  ratio  of  the  variance  of  this  esti- 
U  L  U  L 

mate  to  its  mean  will  be  nearly  that  expected  under  statistical  independence. 

In  the  remainder  of  this  paper,  we  proceed  as  though  the  a(f,  t)  are  independent 
and  distributed  as  ££*,  where  £  is  complex-Gaussian  distributed  with  zero  mean 
and  complex  covariance  matrix  F(f,  t)  (F(f,  t)  is  the  spectral  matrix  at  frequency  f 
and  time  t).  The  primary  consequence  of  this  is  that  sums  of  a(f,  t)  will  be  treated 
as  complex-Wishart  distributed.  Thus,  the  fact  that  a(f,  t)  does  not  have  rank  one 
(except  when  p  =  1)  will  cause  no  apparent  difficulty.  The  use  of  the  complex- 
Wishart  distribution  as  the  distribution  of  a  spectral-matrix  estimate  is  almost 
always  an  approximation.  The  approximation  is  asymptotically  valid  under  various 
circumstances  including  cases  where  the  input  time  series  is  not  Gaussian 
[l2,  p.  220;  32].  There  is  no  asymptotic  result  which  completely  justifies  our 
approximation. 

Our  procedure  for  smoothing  the  elementary  spectral  matrix  estimates  a(f,  t) 
is  a  divisive  clustering  algorithm  that  creates  clusters  connected  in  the  frequency¬ 
time  plane.  Connected  clusters  are  required  because  we  assume  nothing  about  the 
population  spectral  matrix  F(f,  t)  except  smoothness  in  (f,  t).  Thus  there  is  no 
prior  knowledge  implying  that  two  points  in  frequency-time  have  the  same  spectral 
matrix  without  implying  that  the  points  in  between  do  also.  Before  describing  our 
procedure  in  detail,  we  consider  alternatives  to  our  divisive  algorithm  and  our 
dissimilarity  measure. 
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Our  algorithm  starts  with  a  large  rectangular  interval  in  the  frequency-time 
plane  and  partitions  this  into  rectangular  subintervals.  At  each  step,  one  rectan¬ 
gular  subinterval  is  divided  in  two.  There  are  several  considerations  relevant  to 
adopting  this  scheme.  One  advantage  of  a  hierarchical  scheme  is  the  possibility 
that  for  most  inputs  the  gross  behavior  of  each  step  will  be  simple  enough  to  be 
understood  even  though  analysis  of  the  statistical  dependence  of  one  step  on  the 
previous  ones  is  not  tractable.  More  specifically,  our  scheme  which  selects  one 
division  at  each  step  is  simpler  than  a  procedure  that  selects  several  divisions  at 
each  step.  Another  advantage  is  that  hierarchical  schemes  require  fewer  compu¬ 
tations  than  the  most  general  schemes.  The  advantage  of  a  divisive  scheme  is  that 
the  first  steps  are  based  on  large  blocks  of  data  for  which  the  small-sample  biases 
are  less  important.  Since  the  small-sample  biases  of  most  multivariate  criteria 
are  hard  to  evaluate,  these  biases  should  be  avoided  as  long  as  possible.  Another 
aspect  of  this  is  that  for  the  first  step  the  problem  of  statistical  dependence  on  pre¬ 
ceding  steps  does  not  apply  so  that  the  first  step  can  be  analyzed  more  exactly. 
Thus,  the  first  step  can  be  interpreted  for  its  own  significance  such  as  testing  for 
the  uniformity  of  the  population  spectral  matrix  within  the  original  interval.  Some 
alternative  procedures  can  select  partitions  of  the  original  interval  that  our  proce¬ 
dure  cannot  select.  Compared  to  these  alternatives,  our  restricted  range  of  par¬ 
titions  is  an  advantage  because  the  effective  number  of  choices  is  not  too  great  for 
the  available  data  and  because  the  resulting  partition  is  simple  enough  for  visual 
interpretation.  The  disadvantage  is  that  none  of  the  possible  partitions  may  fit  the 
data  very  well. 

In  our  procedure,  we  select  each  division  by  comparing  the  spectral  matrices 
created  from  the  subintervals  on  either  side  of  a  potential  division.  The  criterion 
is  the  largest-root  statistic  adjusted  for  the  variation  of  the  size  of  the  two  sub¬ 
intervals.  At  least  four  ways  to  compare  spectral  matrices  can  be  adapted  from 
tests  for  comparing  covariance  matrices  [  24]  .  For  differences  caused  by  a  rank- 
one  component,  the  largest-root  procedure  is  the  maximum  likelihood  procedure. 
Another  advantage  of  our  choice  is  that  since  spectral  estimates  on  both  sides  of  a 
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potential  division  do  not  have  to  have  full  rank,  a  small  interval  with  a  high  spec¬ 
tral  level  can  be  isolated. 

As  mentioned  above,  one  manifestation  of  bias  in  spectral-matrix  estimation 
is  an  estimate  with  more  significant  principal  components  than  the  population 
spectral  matrix  due  to  smoothing  over  an  interval  that  is  too  large.  This  suggests 
that  intervals  might  be  selected  so  that  the  number  of  significant  components  is 
minimized  in  some  sense.  Algorithms  like  this  have  been  suggested  for  real  mul¬ 
tivariate  data  [9,  p.  289;  13].  These  algorithms  apply  principal-component 
analysis  to  subsets  of  the  data  set  to  discover  whether  the  data  lies  close  to  a 
lower-dimensional  subspace. 

We  complete  the  description  of  our  procedure  by  discussing  the  flow  diagram 
shown  in  Figure  2.  The  inputs  are  a(f,  t)  for  (f,  t)  in  a  large  rectangular  interval. 
This  interval  is  chosen  on  the  basis  of  prior  knowledge  to  be  so  large  that  smooth¬ 
ing  over  a  region  that  extends  from  its  center  to  beyond  its  borders  will  incur 
unacceptable  bias.  In  the  second  block,  the  spectral  matrix  from  each  interval  in 
the  sequence  of  partitions  is  analyzed  using  principal  components  or  techniques 
appropriate  to  the  application.  This  step  does  not  influence  the  sequence  of  parti¬ 
tions,  but  it  does  provide  outputs  helpful  in  choosing  the  right  member  of  the 
sequence.  The  third  block  computes  the  measure  of  dissimilarity  for  all  potential 
divisions.  The  measure  of  dlssimiliarity  for  a  potential  division  is  the  largest  of 
two  one-sided  measures  that  correspond  to  the  two  potential  new  intervals.  We 
formulate  our  procedure  in  terms  of  these  one-sided  measures.  To  specify  these 
one-sided  measures,  let 

s0=  |(f,  t)  lfL<  f  <fw,  tL<  t^tyj  (2.10) 

be  an  interval  already  created  by  the  procedure.  The  new  intervals  that  might  be 

created  by  dividing  are  those  intervals  S  for  which  either  S  or  -  S  is  given  by 

S  =  j  (f,  t)  I  fL<  f  <  ffi,  tL<  t  <  tjj  |  where  fL<  ffi  <  f  or  S  =  j  (f,  t)  I  fL<  f  <  fy, 

t  <  t  <  t  |  where  t  <  t  <  t  .  Let  N  be  the  number  of  points  in  S  ,  let  N  be  the 
L  B  j  L  5  U  0  0 

number  of  points  in  S,  and  let 
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FIGURE  2.  DIAGRAM  OF  THE  PROCEDURE 


.  DIVIDE  MOST  DISSIMILAR  PAIR  OF  SUBINTERVALS 


A  -V 
0  ^  (f,  t)CSf 


a  (f,  t)  ,  A  = 


E 


(f,  t)  C  S 


a  (f,  t). 


The  one-sided  measure  of  dissimilarity  is 


de<S>  = 


if  N<Nq  -  4p 

if  N  >Nq  -  4p  or  Nq<  8p, 


(2.11) 


(2.12) 


where 


d(S) 


1/2 


1 2N  log  N/(Nq^)J  +  2  (Nq  -N)  log  [(NQ-N)/(N0(1-  4*))]  j  if  v  >  N/N 


otherwise 
(2. 13) 


and  4*  is  the  largest  eigenvalue  of  AQ  1  A.  As  discussed  in  the  next  section,  the 

exclusion  of  values  of  d(S)  indicated  in  (2. 12)  is  intended  to  compensate  for  a  small- 

sample  bias.  In  block  4,  the  procedure  finds  the  largest  d  among  those  computed 

e 

from  all  permitted  subintervals  of  all  the  intervals  previously  created.  In  block 

5,  the  subinterval  corresponding  to  the  largest  d  and  the  subinterval  containing  the 

e 

remainder  of  SQ  become  the  new  intervals  in  the  partition.  Equation  (2.12) 
excludes  divisions  of  intervals  with  less  than  8p  points.  Since  SQ  might  be  divided 
into  S  and  -  S  in  two  ways,  d^fS)  is  largest  or  dg(S0  -S)  is  largest,  (2.12)  does 
not  exclude  any  division  of  any  interval  with  NQ  J*  8p. 


Starting  with  the  original  interval,  our  procedure  creates  a  sequence  of 

successively  finer  partitions  until  the  largest  d  is  zero.  There  are  three  guides 

© 

usefbl  in  determining  which  member  of  this  sequence  to  use  in  creating  the  final 
spectral-matrix  estimate.  As  discussed  in  the  next  section,  one  guide  is  the 
sequence  of  values  of  the  maximum  dissimilarity.  In  particular,  regions  for  testing 
the  hypothesis  that  the  first  division  corresponds  to  a  difference  in  the  population 
spectral  matrix  versus  the  hypothesis  that  the  population  spectral  matrix  is  uniform 
are  provided.  Another  guide  is  the  outputs  at  block  2.  These  outputs  can  be  used 
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to  judge  whether  or  not  a  division  reduces  the  number  of  significant  principal  com¬ 
ponents.  Finally,  the  configuration  of  the  intervals  in  each  partition  can  be  used 
along  with  prior  knowledge  to  guide  the  choice  of  partition. 

In  many  applications  of  spectral  analysis,  one  objective  is  the  discovery  of 
narrow  frequency  bands  with  high  spectral  levels.  If  such  a  band  were  in  the 
middle  of  an  interval,  the  set  of  potential  divisions  in  our  procedure  might  not  be 
most  sensitive.  For  example,  comparison  of  the  spectral  matrices  from  the  sub¬ 
intervals  S  =  J(f,  t)  I  f  <f  <  f<f0<fTT,  t<t  j  with  the  spectral  matrices 
from  Sg  -S^  might  be  more  likely  to  isolate  a  narrowband.  Such  comparisons, 
which  would  result  in  two  divisions  instead  of  one,  can  be  included  in  our  procedure. 
Many  of  the  results  in  the  next  section  apply  to  such  an  extension. 


3.  ANALYSIS  OF  THE  CRITERION 


The  core  of  our  procedure  is  the  selection  of  the  two  subintervals  to  which 
correspond  the  most  dissimilar  spectral  matrices.  This  section  is  devoted  to  the 
analysis  of  this  selection  for  a  given  interval.  The  results  are  rigorously  applic¬ 
able  to  the  first  division  of  the  original  interval.  They  can  be  used  for  guidance  in 
interpreting  further  divisions  but  not  with  the  same  rigor  since  the  statistical 
dependence  of  one  step  on  its  predecessors  has  not  been  included.  Because  the 
number  of  alternatives  at  each  step  is  small  compared  to  the  amount  of  data,  the 
guidance  provided  seems  applicable  at  least  for  the  first  few  steps. 

3.1  Maximum  Likelihood 

Our  algorithm  for  selecting  the  most  dissimilar  pair  of  spectral  matrices  is 
the  maximum  likelihood  procedure  for  finding  two  subintervals  whose  spectral 
matrices  differ  by  a  rank -one  component.  Confining  the  set  of  possible  divisions  to 
pairs  of  rectangular  subintervals  as  our  procedure  does,  we  prove  asymptotic  con¬ 
sistency.  Consistency  shows  that  if  the  interval  is  large  enough,  the  bias  does  not 
dominate  the  real  spectral  properties.  Consistency  is  one  indication  that  the  num¬ 
ber  of  choices  is  not  too  large  for  the  amount  of  data. 

Theorem  1:  Let  SQ  be  an  interval  in  the  frequency-time  plane  with  NQ  points, 
where  N„  >  p.  Let  <S  be  a  set  of  subsets  of  with  members  that  are  neither 
empty  nor  have  more  than  N^  -p  points.  Let  the  a(f,  t)  be  independent  and  distri¬ 
buted  as  £|*  where  |  is  complex  Gaussian  with  zero  mean  and  complex  covariance 
matrix  +  uu*  for  (f,  t)  G  S  and  F^  for  (f,  t)  G  -  S.  Let  SG  >S  and  FQ  be  posi¬ 
tive  definite.  The  maximum  likelihood  estimate  of  S  is  the  S  C  that  maximizes 
d(S),  which  is  given  by  (2. 13). 

Proof:  Beside  S,  the  unknown  parameters  to  be  maximized  over  are  the  matrix 
Fq  and  the  vector  u.  The  likelihood  function  is 
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«  ft 


-pNo  r  rN  r  i-(No-N) 

L  =  n  j^det  (F0  +  uu*)  I  I  det  F  J 


j  -  tr  [(FQ  +  uu*)"1  A  +  F’1  (Aq 


-»]!• 


where  A^  and  A  are  given  by  (2.11).  With  probability  one,  the  rank  of  A  is  greater 
than  zero  and  the  rank  of  A^  -  A  is  p.  Maximizing  over  u  with  F^  fixed  gives 

-pN  t  i  “N.  „  |  ^  | 

max^  L  =  II  |det  FQ  j  (a/ N)  exp  j  -  tr  FQ  AQ  +  o  -  N  j , 


where 

a  =  max  j  N,  max  *  ,  w*  F.  A  Fn  w  j .  (3 

j  w*w  =  l  0  0 

Consider  first  the  case  where  A  is  non-singular  and  let  and  4>  be  the  eigen- 

1/2  -1  1/2  J  * 

values  and  eigenvectors  of  A  FQ  A  in  decreasing  order.  Equation  (3.2) 
becomes 


(»/N) 


■exp  iLi=ixj  pj  +  o-N|' 


where 


*  *  .-1/2  .  .-1/2  . 

P.  =  ®.  A  A„  A  9  , 
J  j  0  j’ 


a  =  max  j  N,  Xj  | . 
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Using  the  result 


we  can  maximize  over  \  ,  \  , ...,\  . 

X  z  p 

the  result 


To  maximize  over  4>  ,  6, 

X  A 


we  use 


[n 

j  =  l 


)  I  <  max 
4>,  >  •  •  • »  $ 


[n 


r  JL  -no 

max  I  I  I 

L  * 


4V“,,4p  -j  =  i 


det(A-l/2AA-l/2' 

\  0  , 


P  j  =  1 


N. 


No 

Pj  B(Pj) 


max  p  g(p  ) 
J  *1 


(3.7) 


-N 


0  N 

0  .  . 

«1  gfoy* 


-1/2  -1/2 

where  o^,  »2, are  the  eigenvalues  of  A  '  AQA  '  in  increasing  order. 
The  maximization  of  |  IT  N°  ]  is  a  consequence  of  the  solution  to  the  maximiza¬ 
tion  of  (det  B)N°  exp  (-  tr  B  A  A^  A  [lo].  Thus,  we  obtain  an  expression 
for  the  maximum  of  L  over  and  u  that  is  a  strictly  increasing  fonction  of  d(S). 

If  A  is  singular,  the  proof  must  be  modified  by  initially  replacing  A  by  A+f  I  and 
then  letting  <  go  to  zero  after  maximization. 


In  our  procedure,  the  number  of  possible  divisions  increases  at  least  as  fast 
as  the  square  root  of  the  number  of  data  points.  In  our  interpretation  of  consistency, 
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we  account  for  this  increase.  Thus,  our  theorem  and  proof  require  some  modifica¬ 
tion  of  Wald’s  results  [5,  p.  54],  The  data  sequence  on  which  our  theorem  is 
based  is  composed  of  the  a(f,  t)  fdr  a  sequence  of  nested  rectangles  in  the  frequency' 
time  plane. 


s0(n)=](f.  t) 


w,  t  w<t<t  (n): 
L  U  *  L  U  I 


(3.8) 


We  require  that  as  n-^co,  (f  ^  -  f T^V(tTT^  -  t  ^)- 

U  L  U  L 


(3  where  0  <  p  <  ®  .  The 
parameter  space  that  contains  the  true  parameter' values  and  over  which  the  maxi¬ 


mum  likelihood  estimates  are  formed  is  given  by  F  C  7  .  F  +  uu*  C  7  ,  u*u  >0, 

/„  \  u  v 

and  S  C  S  '  ,  where  7  is  the  set  of  Hermitian  matrices  F  whose  eigenvalues  lie 


between  e  («  >0)  and  ®  and  where  iS 

o« 


(n) 


is  the  set  of  subsets  S  that  along  with 

/yy  V 

S  '  '  is  the  set  of  subsets  incorporated 


-  S  are  non-empty  and  rectangular, 
in  our  procedure.  The  true  parameter  values  are  FQ,  u,  and  S' 
of  generality,  we  let 


,(n) 


Without  loss 


s(n>  ■  j(f. 


t)  f(n)<f<f  <fw,  t(n)<tst(n)i 
’  L  B  U  L  U  I 


(3.9) 


where  as  n  increases  f  remains  fixed  and  neither  f„  -  f,  nor  f,T  -  f  decrease. 

B  ,  .  B  L  U  B 

—  -(n) 

Let  N  denote  the  number  of  points  in  S  , 


Theorem  2:  Let  be  a  neighborhood  of  FQ  and  //^  be  a  neighborhood  of 

Fj^  =  Fq  +  uu*,  where  )I ^  and  are  chosen  so  that  for  all  FQC  //Q  and 

F  G  //  the  largest  eigenvalue  of  F  1  F  is  greater  than  and  bounded  away  from 
a  fa)  a  fa)  /\(n) 

one.  Let  S  ,  F  ,  and  u  be  the  sequence  of  maximum  likelihood  estimates 

over  the  parameter  space  defined  above.  In  the  sequel,  we  drop  the  dependence  on 

n  from  our  symbols  for  the  estimates.  Then,  when  the  a(f,  t)  are  distributed  as 

-(n)  ^ 

above  with  parameters  S  ,  F^,  and  u,  we  have 


Prob  {[s^S^ 

U 

F  l  )L 

U 

F  +  uu*  4 

(l 

0  r  0 

0  r 

1 

(3.10) 
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Proof:  This  is  an  extension  of  ChernofPs  proof  [5].  We  show  that  the  probab¬ 
ility  of  the  event  ^  t  S^j  u  £fq  ^  //Q  J  U  j^FQ  +  uu*  ^  ^  J  happening 

infinitely  often  is  zero.  We  divide  the  sequence  of  occurrences  of  this  event  into 
a  finite  number  of  subsequences  and  then  prove  that  the  probability  of  any  sub¬ 
sequence  occurring  infinitely  often  is  zero.  First,  we  divide  the  event  into  the 

AAA 

cases:  S  contains  the  lowest  frequencies,  S  contains  the  highest  frequencies,  S 
contains  lowest  times,  and  S  contains  the  highest  times.  Each  of  these  cases  must 
be  farther  divided. 


In  order  to  do  this,  we  extend  the  parameter  space  to  a  compact  set.  We  add 
the  points  for  which  u*u  =  0  and  create  the  compact  set  7  *  out  of  7  by  adding  the 
point  at  which  maxw+w_^  w*  Fw  =  ®  .  As  in  [5],  we  create  a  finite  covering  of 
7  *  with  members  //  such  that  for  both  j  «  0  and  j  *  1  if  //  j)  //  then 


<».(//)  =  E  <  inf  log 


j 


'FC)1 


-1  |_.  --1 


(det  F^)  exp  j  -  tr  F^  a(f,  t)j 
(det  F)  1  exp  j  -  tr  F  1  a(f,  t)J 


>0,  (3.11) 


where  the  expectation  is  with  respect  to  the  distribution  for  which  F  is  the  spectral 
matrix.  We  introduce  the  notation 


Q(S,  F,  F)  =  J]  log 


(f,  t)  C  S 


(det  Fy"1  exp  j  -  tr  F."1  a(f,  t)J 
(det  F)  1  exp  J  -  tr  F  1  a(f,  t)j 


Consider  the  events  for  which  S  D  and  F^  +  uu  ^  //^.  For  some  //  ,  we 


(3.12) 


have 


°Slnfs3§(n,  'V'V,/ 

+  Q(sn(S0(h)-S(n)),  F0,  Fj  +Q(S0(n)-S,  F0,  F0)] 


(3.  13) 
(cont) 


>inf 


F^C  n 


Q(S(n),  F, 


1’ 


Fl> 
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+  lnfS  WF  C  7  *  Q(Sn<S0W  •  §<n>)l  V  Fl) 


+  inf  inf  .  Q(S  F  -S.F.F), 
S  F  C  1  *  '  0  0  0; 

0 


By  the  law  of  large  numbers,  we  have  as  n-*-® 


(3.13) 

(cont) 


N  lnfF  cn  Q(S'  \  Fv  F^  —  o^il)  a.s. 


(3.14) 


In  the  second  term  on  the  right  side  of  (3.13),  we  perform  the  minimization  over 
F^  to  obtain 


ln,s  “FjC  7  *  Q<sn<soW  s<n>)l  V  Fi> 


=  infg  log  j  [det  <G/N0) J  G  exp  [  -  tr  (G  -  Ng  I)J  J  , 


(3.15) 


where 


N  1  _  tv  o  -  cn/o  GO  eGOv 


G  =  F0  2-r(f,  t)GS  a(f,  t),  SG  =  Sn(S0'  '  -  S'  '), 


(3.16) 


and  N  is  the  number  of  points  in  S.  We  show  that  N  times  the  right  side  of 
u  G 

(3. 15)  goes  to  zero  almost  surely.  Let  ^  be  the  eigenvalues  of  G  and  let 

1  p 


2  1  (  ) 

-Yj  /2  =  n  j  nq  log  (yiy  -  y  -  nq)  j . 

We  prove  almost  sure  convergence  by  showing  that  [21,  p.  15l] 


(3. 17) 


Pr0b  USUn>k 


£j=i  yj2/2>‘] j 


(3. 18) 


<2  2  Pr°b  2j=iyj/2><  1 ►Oask— >«. 

S  n  ^k  L 


Changing  variables  in  the  density  of  the  eigenvalues  of  G  [17] ,  we  obtain 


21 


Prob  £j=i  yj2/2><]  =FIj  =  i  rt>-J+i)r(NG-j+i) 


pnq  -pn 

ng  . 


/•/ 


Syj  /2>i 


NP  exp  j-N  2j!iyj2/2l  (3’19) 


•riE:i  rij=k+i 

Since  N,  which  appears  in  exponential  In  (3. 19),  goes  to  infinity  at  least  as  fast  as 
the  square  root  of  the  total  number  of  data  points,  (3. 18)  can  be  shown  to  hold. 
Applying  the  same  reasoning  to  the  third  term,  we  see  that  the  right  side  of  (3.13) 
becomes  positive  so  that  the  probability  of  this  subsequence  occurring  infinitely 
often  is  zero. 

The  rest  of  the  proof  consists  of  applying  similar  reasoning  to  the  other  cases. 

In  each  case  we  isolate  a  set  S*  for  which  the  estimated  spectral  matrix  does  not 

belong  to  the  neighborhood  of  the  true  spectral  matrix  and  for  which  the  number  of 

points  goes  to  infinity  fast  enough.  For  this  set  inf  Q  is  shown  to  become  positive 

*  A  ^/n\ 

and  for  the  rest  of  the  interval  inf  Q  is  shown  to  approach  zero.  If  S  DS'  '  and 
Fq+uu*G  then  in  the  set  S*=  £f=fg  +  lj  either  ^  when  !§  =  S^  or 
Fq  +uu*  ^  /?Q  when  is  a  proper  subset  of  S.  In  either  case,  the  previous 
reasoning  applies.  The  case  SCS^  can  be  approached  in  the  same  manner.  In 
the  other  instance  of  frequency  division,  we  consider  first  SDS^  -  S  .  If 
Fq  +uu*  ^  /?Q,  then  the  proof  can  be  based  on  the  set  S*  =  S^n^  -  If 
Fq  +uu*G  then  F  Yl^  or  F^  +uu*^!  Yl^,  so  we  can  use  the  set  S*  =  £f=f^J  • 


The  case  SCS^  -  follows  similarly.  For  time  divisions,  we  treat  the  cases 

S  contains  more  than  half  of  S  ^  separately  from  the  cases  (S  ^  -  S)  contains 
(n) 

more  than  half  of  SQ  . 
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3.2  Distribution  of  the  Criterion 


Consistency  Is  a  property  that  informs  us  about  the  procedure  when  u*u  >0. 
Two  types  of  information  can  be  obtained  by  studying  the  procedure  when  the  popu¬ 
lation  spectral  matrix  is  uniform  throughout  the  interval.  First  is  the  threshold  to 
which  the  maximum  dissimilarity  should  be  compared  in  testing  for  stationarity 
and  uniformity  in  frequency.  Second  is  the  bias  associated  with  the  procedure. 


The  joint  distribution  of  the  largest  and  smallest  eigenvalues  of  *  A  has 
been  derived  [16,  17,  Id] .  From  this  result,  the  joint  density  of  d(S)  and 
d(SQ  -S)  for  a  given  S  can  be  obtained.  The  result  is  given  in  the  following  theorem. 


Theorem  3:  Let  and  X.  be  the  largest  and  smallest  eigenvalues  of  AQ  1  A, 
and  let  N^-p>N>p.  Under  the  null  distribution  (F(f,  t)  =  FQ),  we  have 


Prob  ju<\  <  \  ,<v|  =  c,  det  (b.J, 
P  1  |  1  IJ 


(3. 20) 


where 


i+j-2  N-p  No-N-p 
(x  -  a)  x  (1  -  x)  dx. 


(3.21) 


n  o  r<N0-i  +  l) 

ci~|  li=i  r(NQ-N-i  +  i)  r(N-i+i)  rfr-i+i)  *  (3- 

We  choose  a-  0  or  o=  N/N^  for  different  purposes.  The  distribution  of  for 
N<  p  and  the  distribution  of  \  for  N  >NQ  -  p  can  be  found  similarly  [l7] . 


An  asymptotic  distribution  for  d(S)  can  be  obtained  by  letting  a  =  N/NQ  in 
(3.21)  and  substituting  asymptotic  expansions  for  b_  and  the  gamma  functions. 

Theorem  4:  If  N  and  -N  approach  infinity  in  such  a  way  that  N/N^— ►p  where 
0  <  P  <  1,  we  have  under  the  null  distribution 
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n 


Prob  |o<d(S),  d(SQ -  S)  <  D  j  ~ 

(3.23) 

(n,:,  r.,)-* 

This  distribution  is  shown  in  Figure  3  for  p  =  l,  2,  4,  8,  16  and  32.  Since 

(3. 23)  does  not  depend  on  N  or  N^,  we  see  that  in  a  limited  sense  our  procedure  is 
asymptotically  unbiased  under  the  null  distribution. 

Our  test  for  uniformity  of  the  population  spectral  matrix  within  a  given  interval 
rejects  the  hypothesis  of  uniformity  if  max  r  .  d  (S)  is  greater  than  some  thresh- 

j  © 

old.  Since 

Prob  j  d(S)  or  d(SQ  -  S)  >D  j  <  Prob  j  maxg£-  ^  d(S)  >D  J 

(3.24) 

<  m  Prob  jd(S)  or  d(SQ-S)  >dJ  , 

where  m  is  the  number  of  possible  divisions  of  the  interval,  (3. 23)  gives  approxi¬ 
mate  bounds  on  the  level  of  significance  for  this  test.  As  indicated  below,  N  and 
Nq  must  be  large  for  the  approximation  in  (3. 23)  to  be  close.  Considering  the 
values  achieved  by  the  maximum  dissimilarity  when  the  population  spectral  matrix 
is  not  uniform,  we  see  that  the  upper  bound  given  by  (3. 24)  is  not  too  high  to  be 
useful. 

In  order  to  obtain  a  better  bound  for  the  significance  level  of  our  test  of  uni¬ 
formity  as  well  as  further  insight  into  the  procedure,  we  now  derive  an  approxima¬ 
tion  to  higher-order  joint  distributions  of  the  dissimilarities.  This  result  applies 

only  to  subintervals  S  and  S  -  S  for  which  S  C  S  C  . . .  C  S  .  Also,  the 

K  Ok  mm-1  1 

higher  the  order,  the  more  effort  is  required  in  computing  numerical  values.  How¬ 
ever,  joint  distributions  of  any  order  can  be  used  to  obtain  an  upper  bound  on  the 
significance  level,  with  higher  orders  producing  better  bounds. 

The  following  theorem  gives  the  probability  density  under  the  null  hypothesis 
of  the  eigenvalues  of  A^  *  A^,  k  =  1, . . . , m,  where 
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FIGURE  3.  ASYMPTOTIC  DISTRIBUTION  FOR  VARIOUS  DIMENSIONS 


PROB  (d(S),d(S0-S) * D] 
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VE(f.t)ca  a(f*  *)*  k  _  °» 


m 


(3. 25) 


and  S  C  S  C  •  •  •  C  S  C  S  .  From  this,  the  joint  distribution  of  d(S  ), 
m  m  - 1  1  0  k 

d(S  -St.),  k  =  l,...m  is  obtained  in  Theorem  6.  The  restriction  S  C  S  , 

Ok  m  m  - 1 

C  ...  CSx  means  that  this  theorem  cannot  be  used  to  evaluate  the  dependence  of 
d(S)  for  frequency  divisions  on  d(S)  for  time  divisions.  The  case  of  frequency  and 
time  divisions  appears  to  be  intractable.  The  result  we  derive  applies  only  when 
all  the  roots  of  AQ  1 are  greater  than  all  the  roots  of  Aq  1 A^  + 1  for 
k  =  1, . . . ,  m  - 1.  To  approximate  the  density  for  all  values  of  the  roots,  we  observe 
that  the  probability  that  the  sets  of  roots  overlap  approaches  zero  as  the  number  of 
points  in  +  approaches  infinity. 

Theorem  5:  Let  S  C  S  ,C  . . .  C  S„C  S^;  let  N.  be  the  number  of  points  in 
m  m-1  10  k 

S  ;  and  let  N  ,  N  ,-N  ,  N  -N  . N  "N  ,  N  -N  >  p.  For 

k  m  m-i  m  m-2  m  - 1  2101 

k  =  l,...,m,  let  the  eigenvalues  of  AQ_1  in  descending  order  be  ^iC2,*,*,^kp 
Let  the  a(f,  t)  be  distributed  as  above  with  the  same  population  spectral  matrix 
throughout  SQ.  Then,  when  \  3*  \  ^  +  for  k=  1, ....  m  -  1,  the  joint  density  of 

the  eigenvalues  is 


P^k11*****  =  I  I  ^  , 

11  mp  m  1  li=  1 


N_ “ P  N  -N  -p 

Xml  ‘‘‘Si* 


n»-»  (.,[/.  ,  vWr1 

k=l  jdet[(Xkl  *  X(k  +  1») 


(3.  26) 


where 


nr:iru.  (s,-v<x»i-w- 


c  = 


r(N0-i+i) 


m  1  11=1  I  r(N  -N  -l  +  l)  r(Nm-i  +  l)  T(p-i  +  l) 


1  T  m-l  1 

ilk=1  r<V*W 


(3.27) 
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Proof:  Using  [l7] ,  we  obtain  the  joint  distribution  of 

Bk  =  Ao1/2AkAo1/2>  k=1 . m 


which  is 


p(B  ,  B 

l  l  m 


(3.28) 


(detBm)  m  (del 


■p  <V  (d0t  Bm 


WP 


P  <N  )  Pn(N„-N  ) 
p  m  p  0  1 


(3.  29) 


nr.'. 


N,  -N,  ,  -p 

k  k  +1 


T  (N,  -N.  ) 

p  k  k  +  1 


where  r  (N)  is  the  complex  multivariate  Gamma  function, 


r  (N) 

p 


■^n’.  rpj-i+i). 


(3. 30) 


Using  [17] ,  we  make  the  following  sequence  of  transformations.  For 
k  =  l,  2,...,m,  let  A  =  diag  (\  ,  \  ,...,\  )  and  let  U  be  unitary  matrices. 

K  K1  KZ  Kp  K 

First,  we  transform  B.  to  U,  A,  U,  *  and  then  transform  U„  *B,  U„  to  B.  *  for 

1111  1  k  1  k 

k  =  2, . . . ,  m.  Second,  we  transform  B  '  ’  to  U„  A_  U  *  and  then  transform 

2  2  2  2 

U  * B.  ^ U.  to  B. ^  for  k  =  3, . . . , m.  We  continue  this  until  we  transform  B  ^ 
2  K  2  K  m 

to  U  A  U  *.  Integrating  over  U.  for  the  case  \,  ^  \ .  we  obtain 

m  m  m  k  kp  (k+l)l 


/  \Nm"P/  \N0"Nl'P 

P(Hl . W=Cm(detAm)  r,(I-AX>) 


r  /  .  n  -  n. 


n.;. 


(3.31) 

(cont) 


(dU) 


rii=i  [r<vNk+i-,+i>i'*-i+i>/r<v,w] 
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n:.n::;ru. 


(Xki"Xkj)  ’ 


(3.31) 

(cont) 


where  (dU)  Is  the  Invariant  measure  on  the  unitary  group  U(p)  normalized  to  make 
the  total  measure  unity  [14] .  Note  that  when  ^kp<  +1^  the  range  of  integra¬ 
tion  of  U  is  limited  by  the  requirement  that  A  -  UA  U*  be  positive  semi- 

k  k  +  1 

definite. 

Using  [14] ,  we  can  evaluate  the  integrals  over  U.  We  have 


in“p 

det  (I  -  AUBU*)  J  (dU)  =  1 FQ  (-n  +p;  A,  B) , 


(3. 32) 


where  A  -  diag  (0^, . . . ,  o^),  B  =  diag  (p^, . . . ,  p^),  and  is  a  hypergeometric 
function  of  matrix  argument  defined  in  [14]  .  Using  the  definitions  provided  in 
[14]  ,  we  obtain 


1FQ(-n+p;  A,  B)  = 


£k  =  0  ^k  >  ...^k 
1  P 


n.p. 


r(p-i  +  l)  (-n+p-i  +  l). 


r(k.+p-i  +  l) 


det 


/  kj+P-Jv  ,  Mp-j, 

K  )  det  \^i  ) 


det  |a.P  ^jdet^p^ 

where  the  sum  over  k^>  . . .  ^  k^  is  the  sum  over  the  partitions  of  k  and 
(n)^  =  (n)  •  (n  +  1)  ...  (n  +  j-1).  To  sum  (3.33),  we  need  the  identity 

/  kj+p'j\  /  kj+p-j\ 

det  ^  jdet^pj  J 


(3. 33) 


*  J  f  •  •  •  f  J 


sign  (r  , . . . ,  r 


>*  ^si . ®p  FI i=  1  (°r(s.)  ^sj 


k  +p-i 
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(3.34) 


where  the  summations  are  over  the  permutations  of  (1,  2,... ,p),  r(s)  denotes  the 

result  of  successive  permutations  of  i  by  r  and  s,  and  sign  (r  , . . .  ,r  )  is  +1  (-1) 

*  P 

if  the  permutation  is  even  (odd).  By  summing  the  binomial  series,  we  obtain 


det  ((a  jP  "  J )  det  (P  "  ^ )  j  +  P5  A,  B) 


[n 

Li  =  1 


r(p  -  i  +  i) 


^  sign  (r1#  •  •  • » rp) 


ri’-**’rp 


p  r  (-n  +  p  -  i  +  1) 


iC  2  £  fl  r^  +  p-l  +  1)  (°r(s.)  pSl) 


k=0  k  J* . . .  >k  s  i  ■ .  ■  i  s  i=l 

1  pip 


k.  +  p  -  i 


(3. 35) 


n 

i  =  1 


r(p  -  i  +  i) 


(-n  +  1) 


p  -  i  J 

P 


ft 

2  8lgn(ri . V  n  (‘‘"r  pl)”"1 


r  > . . . ,  r 
1  P 


i  =  1 


=  (-1)P(P  "  1)72  Y\  [r(P  -  i  +  1)  Hn  -  i  +  l)/r(n)J  det  ((1  -  a^)*1  1 J  . 

i  =  1 

Substituting  (3.32)  and  (3.35)  into  (3.31)  completes  the  proof. 

Theorem  5  provides  a  formula  for  the  density  that  is  exact  when  the  roots  from 
different  intervals  do  not  overlap.  We  now  integrate  this  density. 


Theorem  6:  If  the  conditions  of  Theorem  5  hold  and  if  v  >  u  >  v  >  u0. .  .> 

1  1  Z 


v  >  u  ,  then 
m  m 


Prob  j  vk>XRi . \kp>  i^;k  1 . mj  cm  det  (a{j) 


(3.  36) 
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where 


,  -T1  f2  ...T”  (X  -»/-l(x  -./-j 

11  >1  >2  Jum  m  1  12 


N  -p 

Jx  m 

m 


(3. 37) 


/  \N0’Ni'P  rTm-1  /  \  \"k  A'k  +  1  * 

*(1"Xl)  1  |k=l  (k  k+l)  dXm***dXl* 

The  constants  and  can  be  chosen  at  our  convenience. 

Proof:  Since,  for  any  i,  j,  and  k,  X^  and  X.  can  be  interchanged  without 

changing  the  value  of  the  expression  for  p(X,, , ....  X  )  given  in  (3. 26),  we  have 

li  mp 


N  -  N.  - 1 


fV  1  /*11  Am(p-1) 

/u  i  J^'J^  11 . 


X  )  d X  ...  dX 

mp  mp  11 


■  h-rn?  -l': . 


(3.38) 


X  )  dX  . . .  dX  . 
mp  mp  11 


Using  the  identity 


rif=i  ru,  <^i-w=det  ((*tai-"i)p'j)’  <3-39) 


and  expanding  the  determinants,  we  have 


n,=P: 


Xmi  det((X(m-l)i“\nj 


N  -N  -1\ 
j  m  - 1  m  j 


•nr::1  n ,.r+,  <^,-w 


(3. 40) 


=  Er . r  Es . s  sl«”<8i . V 


1  P  1 


np  i 

1  li=l  mrl 


N  -p 
m 


(X(m  -  l)i  Xmri  ) 


N  -N  -1 

m  - 1  m  .p-s, 

Xmrl~a1  l* 
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where  r ...... r  and  s„ , . . . ,  s  are  permutations  of  1 . .  Integrating  (3. 40) 

1  pi  p 

over  u  <  X.  v  ,  i  =  l,...,p,  we  obtain 
m  mi  m 

r(P+1)2s1 .  slgn(8l . vn.=  i[/umm  xm  "  P 


N  _ “N  -1 
m 


(w*.)"-'1  "m 


p-s, 


dX. 


ml 


(3.41) 


v  N  -p 
.  — .  _  .  ■  ■  m  m 

=  r(p+l)  det  I  I  x. 

1  fu  m 
m 


(X(m-l)i~\n 


N  -N  -1 
m  - 1  m 


(v-.r 


d\ 


m 


Repeating  this  process  proves  the  theorem. 


An  approximation  for  the  distribution  of  \  ,  k  =  1, . . . ,  m  valid  for  large 

N  -N  ,  N  -N . N  ,  - N  ,  N  can  be  obtained  by  extending  the  lower  limits 

0112  m-lmm 

in  (3. 37)  into  the  region  where  the  density  is  not  valid.  When  the  contribution  from 
this  region  is  small,  this  approximation  is  useful.  We  obtain 


Prob  |vk>>^kl;  k  =  l,...,m  j  ~  cm  det  (a'  ), 


(3.42) 


where 


vf-c-  '/;■  K-.r' (‘.-.r 


(3.43) 


*  \ 


N  -p 
m 


N  -  N,  - 1 


m 


•  d  . . .  d  X.  . 
1  m 


It  can  be  shown  that  (3. 42)  gives  the  correct  total  probability,  Prob  1  >  X.  ; 

,  I  K1 

k-l,...,m,  -1. 
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'  0 


For  the  case  p  =  l,  (3. 36)  and  (3.37)  can  be  recognized  as  the  probability  that 

the  sample  quantiles  of  order  N„/N  .  N_/N,. , . . . ,  N  /N„  lie  within  the  given 

1  0  20  mO 

limits  for  independent  samples  from  a  uniform  (0,  1)  distribution.  Thus,  the  rela¬ 
tion  of  our  procedure  to  the  test  presented  in  [8]  for  a  single  time  series  can  be 
seen.  This  test  is  based  on  the  maximum  of  |  -  k/N^|  over  1<  k<  NQ.  Thus, 

f  <  f  <f  +k  -  1,  t=t  {  and  uses  a  different  function  of 
L  L  L ) 


this  test  chooses  =  j  (f,  t) 

as  the  measure  of  dissimilarity. 


Setting  a,  =N  /Na  and  a  =N,/N.  in  (3.37),  we  can  obtain  an  asymptotic 
l  m  0  2  1  0 

expansion  for  a^  and  thus  an  asymptotic  distribution  for  d(S^),  d(SQ-S^), 
k  =  1, . . . ,  m.  The  expansion  of  a„  is  nearly  the  same  as  that  obtained  in  the 
asymptotic  distribution  of  sample  quantiles  [6,  p.  201] . 


In  order  to  obtain  moments,  we  Consider  an  equivalent  formulation  of  our  pro¬ 


cedure.  Let  A.'  ,  X'  , . . . ,  X*  be  the  eigenvalues  of  1 A 
kl  k2  kp  Ok 

order  and  let 


taken  in  a  random 


ykj  = 


(WNo) 

Lkiog 

\  i 

|x'  -n,/n  I 

NX! 

1  kj  k  0| 

[ 

L  0  kj  J 

+  2(NQ  -  Nfc)  log 


VNk 

Wj 


-.il/2 


(3.44) 


Our  procedure  is  equivalent  to  selecting  S  G  S  if 

tv 


ykjl  =  m“k.j  4) 


value  of  j.  We  obtain  expressions  for  the  second  moments  of  y,  . 

kj 


for  some 


Theorem  7:  Let  y^,  y ■  be  defined  as  in  (3.44).  If  NQ,  N^,  N2>  approach 
infinity  in  such  a  way  that  N  /N  — ►  p  and  N  /N  — ►  P  ,  0<  p  <  pel,  then  we 

1  v  1  «  U  «  6  1 

have 
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Eyy  E52j~° 


„  2  _  2 
EVEy2j  ~P 

Eyliyij'  Ey2ly2J~-1  i*i 


(3.  45) 


Ey!iy2j  ~  P_1  [p2 (1- 


1/2 


Proof:  We  first  derive  the  moments  of  X.^  and  then  use  the  expansion 

ykj'  [v«3k<i-pk»]1/2  <w- 

Proceeding  as  we  did  in  the  proof  of  Theorem  6,  we  obtain 


(3.46) 


P(\'u....,  \’2p) 


cJr(PU)l-2(-l)IJ«,-1)/2V  Q  T  V 

*  L  J  ql,***,qp  rl . r 


N  -i 


sign(si . v  n,:  w 1  (vs) 


W1 


(3.47) 


/  VN  -N, -».  ) 

j 


We  integrate  applying  the  approximation  used  in  (3. 42)  and  (3.43)  to  obtain 


E^„'-'-l>P<P'1,/2c2[r(P  +  l)] 


-1 


det 


qi,,,,,qp 


T(N2+6(n,  q^-i+l)  r  (N^  -  Ng)  r  (Nq  -  -  j  + 1) 

r(NQ+  6(n,  qt)-i-j  +  2) 


(3.48) 


W^ninp-.VjE 


(N  -k  +  1)  •  det 

Ct 


i  =  1 


k  =  1 


r(N0-j  +  i) 


r(NQ+  6(i,  k)  - 1  -  j  +  2) 
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where  6(i,  j)  is  the  Kronecker  delta.  Making  use  of  the  Identity 
Nq  det [r(NQ - i + 1)/ r(NQ  +  6 (i,  l)-i-j+2)]  = 


(3.49) 


=  det  [(N0+6(i,  ll-i)3”1]  =  (-1)P(P“1)/2  p  r(p-i  +  l), 

we  see  that  E  X.'  ~  N  /N  .  The  derivation  of  E(1  -  X*  .),  E(  Xr  X*  ),  and 
2j  2  0  lj  21  2j 

Etfi-X^)  (1-  X'  ))  proceeds  similarly. 

Consider  finally  the  derivation  of  EfX'^O--  X'  )),  where  without  loss  of 
generality  we  let  i  =  j  =  l.  Integrating  and  replacing  r^  by  rfs^),  we  obtain 

4,21(1-Kil))  =  (-DP<P'1)/2  =2  [r(P+1>]"2 


*  Y!  „  YL  „  V.  o  sign  (s  ,...,s  ) 

i  p  i  p  i  p 


n 

i=  1 


r<N2+6(l,  q^-i+D  r^N1  “ Ng)  r (NQ - Nx  +  6 (1,  r(si))-si  +  l) 
r(NQ+  6(1,  qt ) +  6 (1,  r(s.))  -  i  -  s.  +2) 


(3.50) 


=  (-l) 


p(p-l)/2 


-1 


kh "  t  Ek£iZk-?i 

p 


r(NQ-i  +  i) 

r(N0+6(i,  k)+  6(j,  kf) -  i- j  +2) 

In  (3. 50),  the  determinant  is  non-zero  only  when  k  =kf  =  1.  The  determinant  in 
(3.50)  equals  *  times  the  following  determinant  which  we  evaluate  by  expanding 
about  the  elements  of  the  first  column  and  evaluating  the  minors  using  results 
similar  to  (3. 49).  We  obtain 


(N2-k+l)  (NQ  -Nx  -k*  +l)det 
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det 


r(N  +6(if  1)-1+1)  ) 

r(NQ+  6(i,  1)+  6(j,  l)-i-J+2)j  =  Jl1^1  (N0+  6(1’  1)_1) 


m  =  1 


H) 


m  - 1 


(NQ+  6(m,  l)-m  +  l)  (Nq+  6(m,  l)-m) 


n  f- 11  n 

i  f  m  j  ^  m 


+  D 


=  (-l)P(P"1)/2  n[i  =  l  [r(P-i  +  l)  (NQ- 6(i,  l)-i)J 


(3.51) 


E_  V  P  (-l)P"m/P_1\ 

r(p)  Z-fm  =  l  v  ’  I  m-1  J 


(N0“m  +  1)  (N0-m) 


7Fi)Sm=P0  H>P""(P)  |(N0-m.l)(N0-m) 


-1 


-1 


The  summations  on  the  right  side  of  (3. 51)  are  evaluated  using  the  identity 


EP  /  TvP-m 
m  =  l  <"1) 


(-•’>)  i 


(NQ-m+l)  (NQ-m) 


-1 


rip  Vp-1,,  ,p " i ,  , 

II  x  (1-x)  dxdy 

yo  Jo 


(3.52) 


=  r(p)  r(N0-P+i)/r(N0) 


Ep 

m  =  l 


(Nq ~P  +  m)  (Nq “P  +m - 1) 


-1 


and  the  equivalent  identity  obtained  by  replacing  p  by  p  + 1.  Equation  (3. 45)  follows 
from  the  0(1)  and  0(Nq  *)  terms  of  the  result. 


3.3  Small  -  Sample  Bias 

We  would  like  to  have  a  procedure  that  under  the  null  hypothesis  choses  any 
division  with  equal  probability.  Alternatively,  a  procedure  that  tends  to  divide  the 
interval  into  equal-sized  subintervals  is  preferable  to  other  tendencies.  Neither  of 
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these  properties  hold  for  our  procedure,  and  we  have  excluded  certain  values  of 
d(S)  from  our  selection  to  improve  its  behavior.  As  shown  by  the  form  of  the 
moments  in  (3.45),  our  procedure  does  not  select  divisions  with  equal  probability 
even  asymptotically.  There  is  also  a  small-sample  bias.  We  investigate  this 
problem  by  studying  the  univariate  and  bivariate  distributions  derived  above  and 
by  simulation. 


Under  the  null  hypothesis,  the  maximum  of  d(S)  tends  to  occur  for  large  N.  In 
the  extreme,  we  see  that  d(S)  =  co  for  N  >N^-p.  One  indication  of  this  bias  is 
given  in  the  following: 


Theorem  8:  When  N  >Nq/2,  we  have  under  the  null  hypothesis 
Prob  Jd(S)^d(SQ-S)>oj  >  Prob  jd(SQ-S)^d(S)>  oj  , 


(3.53) 


where  a  is  any  constant. 


-1 


(1) 


Proof:  Let  X  ,  X  , . . . ,  X  be  the  eigenvalues  of  A  A  in  decreasing  order. 

J.  Z  p  0 

Dropping  the  unneeded  subscript,  let  y^  be  defined  as  in  (3. 44).  Since  this  function 
is  monotonic,  we  can  consider  X  ^  to  be  a  function  of  y^,  =  X(y^).  Let  j 

X^Vj)  *  dXj/dy^.  Using  [l7]  ,  we  obtain 

Prob  |d(S)>d  (SQ-S)>  o|  =  Prob  |  -yp  >  aj 


f 


PN(yi)dyi' 


(3.54) 


where 


pN(yi)  r(p) 


IP 


(t)  (V) 
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J-h  'hi  "p  ‘TI-JMXj 
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-yx2/2 


(i) 


(3.  55) 
(cont) 
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•  rii=2  n5)  dy2---  %• 


(3.55) 

(cont) 


If  in  (3. 44),  we  exchange  N  with  NQ  -  N,  we  obtain  another  function  of  which  we 
denote  4  =  4*(y  ).  Since 

J  J 


\(y^)  =  1  .  4-(-y^);  \(1)  (y^)  =  4^  (-y^) 


(3.  56) 


(1) 


we  can  obtain  an  expression  for  Pnq  _n  6^)  by  replacing  X^  with  4^  and  X^  with 
4-^  in  (3. 55).  We  complete  the  proof  by  showing  that  pN  (y^)  >  pNq  _  N  (y^).  The 
curve  (X(y^),  4*  (y^))  nasses  through  (N/N^,  (Nq  -  N)/N^  and  (1,  1),  and  at 


y^=0,  X  a)  =  +  a>.  Since 


4,  (1) 

'  Vi-N 

V-v 

Tv  ^ 

i 

Vr<VN> 

x1a-^1) 

_  «N 

(3.57) 


we  see  that  >1  if  N()X1  -N  >1*^  -  (NQ  -  N)  and  4  ^/X^^l  if 


^i>xi‘  Thus»fory1^°* 


(3.58) 


This  completes  the  proof. 

Figure  4  gives  Prob  |  d(S)<  D  j  and  Prob  |  d(Sj)  and  d(S^)<D  |  under  the  null 

hypothesis  for  N  =128  and  various  N,  N  ,  and  N  .  The  horizontal  scale  is  a  non- 
U  1  z 

linear  function  of  D  chosen  so  that  asymptotic  distribution  of  d(S)  plots  as  the 
identity  function.  We  make  three  observations.  First,  the  probability  of  d(S) 
exceeding  a  given  value  increases  with  N  as  does  the  probability  of  d(S  )  or  d(S  ) 

X  z 

exceeding  a  given  value.  This  fact  agrees  with  the  simulation  result  that  the 

largest  d(S)  tends  to  occur  for  intervals  with  large  N.  Second,  for  large  D,  the 

sum  of  the  probability  of  d(S  )  >D  and  the  probability  of  d(S  )  >  D  is  not  much  larger 

X  z 
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than  the  probability  of  d(S  )  or  d(S  )  >D.  This  indicates  that  the  univariate  distri- 
buttons  can  be  used  to  obtain  an  adequately  tight  upper  bound  on  the  level  for  testing 
uniformity.  Finally,  for  a  given  D  the  value  of  the  distribution  for  N  =  64  is  smaller 
than  the  value  of  the  asymptotic  distribution.  This  fact  and  the  results  for  the  other 
values  of  N  indicate  that  the  asymptotic  distribution  does  not  give  true  upper  bounds 
on  the  level  of  significance.  This  also  indicates  that  our  procedure  might  tend  to 
divide  small  intervals  in  preference  to  large  intervals. 

To  obtain  results  on  the  small-sample  bias  that  include  all  effects,  we  per¬ 
formed  500  replications  of  one  step  of  our  procedure.  The  data  were  identically- 
distributed,  8-dimensional,  complex- Gauss lan  vectors  for  the  interval 
Sq  =  j  (f,  t)  |  l<f<16,  l<t<8  j.  Figure  5  shows  the  histogram  of  divisions  chosen. 
For  the  time  division  t-j,  the  shaded  contribution  to  the  relative  frequency  resulted 

from  d  (S)  being  largest  forS=  ]  1  <  f  <16,  1<  t  <  j ! ,  and  the  unshaded  contribution 

6  I  I 

resulted  from  de(S)  being  largest  for  S  =  j  1<  f <16,  j<  t<8  j .  The  results  for  the 

frequency  divisions  are  displayed  similarly.  The  tendency  for  the  maximum  of 

d  (S)  to  correspond  to  an  interval  with  large  N  is  clearly  seen, 
e 

Our  procedure  excludes  d(S)  when  N  >NQ  -4p.  This  is  why  divisions  t-7,  f-14, 
f-15,  and  f-16  have  no  shaded  contribution  and  why  division  t-1,  f-1,  f-2,  and  f-3 
have  no  unshaded  contributions.  Without  this  exclusion,  the  simulation  results 
were  much  more  biased.  The  two  intervals  with  N  =  120  were  chosen  473  times  out 
of  500. 

The  simulation  results  can  be  used  to  illustrate  an  upper  bound  on  the  level  of 
significance  for  our  test  of  uniformity.  In  the  simulation,  9%  of  the  500  values  of 
max  c  d  (S)  exceeded  6.5.  Using  (3. 20),  we  obtain  0.12  for  the  sum  of  the 
probabilities  of  [d(S)>6.5]  for  the  36  intervals  examined  in  the  simulation.  Thus, 
this  upper  bound  is  apparently  tight  enough  for  most  purposes.  If  we  multiply  the 
asymptotic  probability  of  [d(S)  >6. 5]  by  22,  the  total  number  of  divisions,  we 
obtain  0.03.  Thus,  in  this  case,  the  asymptotic  distribution  is  for  enough  above 
the  true  distribution  for  most  N  that  it  underestimates  the  significance  level. 
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FIGURE  5.  SIMULATION  FOR  NULL  DISTRIBUTION 
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DIVISION  SELECTED 


The  usefulness  of  an  upper  bound  based  on  (3. 20)  can  be  expected  to  vary  with  p. 


In  particular,  since  (3.45)  seems  to  indicate  that  the  dependence  of  d(S  )  on  d(S  ) 

J.  dt 

decreases  with  increasing  p,  such  bounds  may  not  be  as  tight  for  smaller  p. 
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4.  SIMULATION 


Beside  the  analytical  complexity,  the  broad  range  of  alternatives  to  the  null 
hypothesis  makes  precise  specification  of  the  performance  of  our  procedure  in  all 
situations  nearly  impossible.  In  this  section,  we  present  the  results  of  a  simula¬ 
tion  for  a  simple  situation  for  which  the  correct  behavior  is  easily  seen.  A  simu¬ 
lation  of  a  realistic  sonar  situation  has  been  presented  [20].  The  simulation  pre¬ 
sented  in  this  section  shows  that  the  sensitivity  of  our  procedure  increases  with  p, 
the  number  of  time  series  analyzed.  This  suggests  that  determination  of  smooth¬ 
ing  for  each  single  time  series  in  a  multiple  time  series  followed  by  application  of 
the  result  to  the  analysis  of  the  multiple  time  series  is  not  as  sensitive  as  our  pro¬ 
cedure.  This  simulation  also  provides  information  on  the  magnitude  of  the  differ¬ 
ences  our  procedure  will  detect. 


The  data  we  created  were  16-component,  complex-Gaussian  vectors  for  the 
interval  SQ  =  j  (f,  t)  I  1  <  f  <  40,  1  <  t  <  16  J  .  The  population  spectral  matrix  was 

u  u  *  +  uu„  *  +  uu„  *  +  u  u  *  +  I  for  1  <  f  <  8 
11  22  33  44 

u„u  *  +  u„  u  *  +  u  u,  *  +  I  for  9  <  f  <  16 
2  2  3  3  4  4 


F  = 


u  u  *  +  u  u  *  +  I  for  17  <  f  <  24 
3  3  4  4 

u  u  *  +  I  for  25  <  f  <  32 
4  4 


(4.1) 


I  for  33  <  f  <  40, 


where 


\ 

/ 

(4.2) 
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Because  u  ,  u  ,  u  ,  and  u  are  mutually  orthogonal,  the  population  eigenvalues 

1  a  o  4 

are  easily  determined. 

We  created  20  independent  sets  of  data.  We  applied  the  procedure  to  each 
single  time  series,  to  consecutive  pairs  of  time  series,  etc.  to  obtain  320  inde¬ 
pendent  relications  for  p  =  1,  160  independent  replications  for  p  =  2,  80  for  p  =  4, 

40  for  p  =  8,  and  20  for  p  =  16.  The  results  for  different  values  of  p  are  dependent, 
a  fact  which  does  not  make  this  simulation  unsuitable  for  determining  the  increase 
in  the  sensitivity  with  p. 

The  results  of  the  simulation  are  shown  in  Figure  6.  We  denote 
!  (f,  t)  I  f  <  f  <  f  ,  1  <  t  <  16 }  by  [f  -  f  ].  Ideally,  our  procedure  should  first 
separate  [l-8]  from  [9-40],  then  separate  [9-16]  from  [l7-40],  then  separate 
[17-24]  from  [25-40]  and  finally  separate  [25-32]  from  [33-40]  .  Each  of 
these  divisions  should  occur  with  the  lower-frequency  side  having  the  extra  com¬ 
ponent.  This  ideal  behavior  occurred  four  times  when  p  =  16. 

In  Figure  6a,  the  relative  frequencies  for  the  first  step  are  shown  for  p  =  1 
and  2.  In  the  first  step,  the  maximum  of  de(S)  always  occurred  for  an  S  of  the  form 
|  (f,  t)|l<f  <f^  1  <  t  <  16  j .  The  possibilities  division  in  time  and  the  maxi¬ 
mum  occurring  for  the  high-frequency  side  did  not  occur.  For  p  =  2,  the  correct 
division  was  made  in  almost  half  the  replications.  For  p  =  1,  the  divisions  were 
more  spread  out.  For  p  =  4,  the  correct  first  division  was  always  made  with  one 
exception. 

In  Figure  6b,  we  show  for  the  second  step  the  sum  of  the  relative  frequencies 
of  the  interval  pairs  J  [l-7  ] ,  [8-^]J,  |[l-8],  [9-f]|,and  j[l-9],  [l0-^]J. 
For  p  =  4,  all  replications  had  this  form,  one  having  been  created  by  first  select¬ 
ing  [l-16]  and  then  dividing  it  into  [l-8]  and  [9-16],  Note  that  when  p  =  4, 

[7,  8  or  9-16]  occurred  in  over  half  the  replications.  For  p  =  2,  131  of  the  repli¬ 
cations  had  the  form  appropriate  to  Figure  6b,  29  of  these  having  been  created  by 
first  selecting  [l-'f']  and  then  dividing  it. 

In  Figure  6c,  we  show  for  the  third  step  the  sum  of  the  relative  frequencies  of 

the  interval  triples  j [l-8],  [9-15],  [l6-f]j,  j [l-8],  [9-16],  [l7-^]J,  and 
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FIGURE  6.  SIMULATION  FOR  FIVE  SUBINTERVALS 
WITH  DIFFERING  SPECTRAL  MATRICES 


INTERVAL  SELECTED  -  f 
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{  [l_8] ,  [9-17] ,  [l8-^]  | .  For  p  =  4,  the  results  were  not  even  approximately 
correct  more  than  half  the  time.  For  p  =  8  and  16,  the  first  division  always 
separated  [l-8]  from  [9-40].  For  p  =  8,  one  replication  did  not  have  the  form 
appropriate  to  Figure  6c.  For  p  =  16,  they  all  did.  Note  that  24  occurred 
over  half  the  time  for  both  p  =  8  and  16.  Even  for  p  =  16,  the  results  of  the  fourth 
step  were  usually  incorrect. 

A  few  generalizations  about  the  performance  of  our  procedure  in  the  situation 
simulated  seem  noteworthy.  We  see  that  doubling  p  does  not  quite  double  the  sen¬ 
sitivity  and  that  when  u*  FQ  1  u  (which  equals  u*u  when  p  3*  4)  is  greater  than  one, 
the  procedure  performs  reasonably.  The  results  for  p  -  16  may  indicate  that 
increasing  p  may  not  increase  sensitivity  when  the  size  of  frequency-time  interval 
is  too  small. 

The  computer  time  needed  for  the  simulation  is  of  interest.  Each  replication 
consisted  of  five  cycles  through  the  procedure  (creating  six  intervals).  On  a  CDC 
6700,  each  replication  took  on  the  average  0. 60  seconds  for  p  =  1,  1. 73  seconds 
for  p  =  2,  4. 52  seconds  for  p  =  4,  15. 8  seconds  for  p  =  8,  and  71. 7  seconds  for 

p  =  16. 
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